\documentclass{article}
\usepackage{xltxtra}
\usepackage{ctex}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{geometry}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{float}
\usepackage{geometry}
\usepackage{booktabs}
\usepackage{tikz}
\usepackage{listings}
\usetikzlibrary{shapes.geometric}
\usetikzlibrary{shapes.arrows}
\usetikzlibrary{positioning}
\usetikzlibrary{graphs}
\usetikzlibrary{quotes}
\usetikzlibrary{arrows}

\makeatletter
\newcommand{\Rmnum}[1]{\expandafter\@slowromancap\romannumeral #1@}

\makeatother
\geometry{a4paper,top=1.5cm,bottom=1.5cm}
\title{Report for project.}
\author{张皓祥 \\ 3200102536 强基数学2001}
\date{}

\begin{document}
\maketitle

% 构建一个能够实现构造牛顿插值多项式并做运算的项目
\section{项目构建思路}

本项目旨在实现节点插值拟合曲线，分为了根据多项式插值拟合 (ppForm) 和 B 样条插值 (B-splines) 两种实现方式。

其中，ppForm 的插值方式实现了 $S_1^0$ 的线性插值， $S_3^2$ 中三次函数插值的 
Complete, Specified, Nature, Not a knot, Periodic 
五种边值条件。而 Bspline 的插值方式实现了 $S_1^0$ 和 $S_3^2$ 中的 Complete, Specified, Nature, Periodic 四种边值条件，
同时根据题目要求，实现了基于 Complete 条件的 Cardinal 类型的 $S_3^2,S_2^1$ 类型的插值计算。

实现两种插值计算方式分别依赖于两个求解类 \verb|class ppForm| 和 \verb|class Bspline|. 
其中，为了构建 ppForm 插值，规范了仿函数类 \verb|class Function|, 同时由于该种插值方式基于多项式拟合，建构了多项式类 
\verb|class Polynomial| 用于规范多项式运算。

而为了构建 Bspline 插值，除了规范仿函数类，更重要的定义出 B 样条基函数。
为此，该项目首先构建了单个 B 样条基函数类 \verb|class Bbasis|, 然后运用单个基函数，构造了一个能够根据给定插值节点
和次数构造整体 B 样条基函数的类 \verb|class Bbase|. 

编写时，两个求解器公用的结构和宏定义被写入
一个库 \verb|InterpolationStruct.h| 中。整个项目整体构建思路如下：
\bigskip

\begin{tikzpicture} [node distance = 20pt, bend angle=45, auto]
    \centering
    \node [rounded corners, draw]                      (fun)    {class Function};
    \node [rounded corners, draw, left=45pt of fun]    (pol)    {class Polynomial};
    \node [rounded corners, draw, right=45pt of fun]   (basis)  {class Bbasis};
    \node [rounded corners, draw, below=20pt of basis] (base)   {class Bbase};
    \node [rounded corners, draw, below=45pt of pol]   (pp)     {class ppForm};
    \node [rounded corners, draw, below=20pt of base]  (bspline){class Bspline};
    \node [rounded corners, draw, below=20pt of fun]   (inter)  {InterpolationStruct.h};
    
    \draw [->] (fun.west)    to (pol.east);
    \draw [->] (fun.east)    to (basis.west);
    \draw [->] (pol.south)   to (pp.north);
    \draw [->] (basis.south) to (base.north);
    \draw [->] (base.south)  to (bspline.north);
    \draw [->] (inter.west)  to (pp.east);
    \draw [->] (inter.east)  to (bspline.west);
\end{tikzpicture}

% class Function
\subsection{class Function 仿函数抽象类}
项目中需要涉及到大量函数的运算，除了多项式，B 样条基函数本身是一类函数，在做插值过程中也会碰到各种函数。为此，本项目中定义了
一个仿函数抽象类。其一方面可以规范所有相关函数的公有特性，另一方面可以为所有仿函数类的传递提供一个统一媒介。
该仿函数类存放于文件 \verb|Function.h| 内，其中仅定义了三个函数：
\begin{verbatim}
    virtual double operator() (const double & _x) const = 0;
    virtual double diff(const double & _x) const;
    virtual double second_diff(const double & _x) const;
\end{verbatim}

其中，重载操作符 () 的操作定义为一个纯虚函数，用于规范所有仿函数类都应该对 () 进行重载，并且重载规定输入和输出
样式，实现统一的调用接口。
而函数 \verb|double diff| 和 \verb|double second_diff| 则是用于对函数进行求一阶导和二阶导。
该函数作为虚函数可以被覆盖继承。

% Interpolation.h
\subsection{公共库 InterpolationStruct.h}

该库函数旨在定义求解器需要的共同结构，以及提供一些必要的宏定义。

\noindent $\bullet$ 调用所有需要的辅助库

调用了两个求解器都需要使用的库函数，包括 \verb|<iostream>|, \verb|<vector>|, \verb|<algorithm>|, 
\verb|<cmath>|, \verb|"Function.h"| 以及用于解线性方程组的 \verb|Eigen| 库。

\noindent $\bullet$ 定义插值类型 
\lstset{language=C}
\begin{lstlisting}
# define S01 1    // 线性插值
# define S32 3    // 三次插值
# define S21Car 2 // Cardinal S21 二次插值
# define S32Car 4 // Cardinal S32 三次插值
\end{lstlisting}

\noindent $\bullet$ 定义边值条件: \verb|Complete|, \verb|Specified|, \verb|Nature|, \verb|NotAKnot|, \verb|Periodic|. 

\noindent $\bullet$ 定义节点单元: \verb|struct point|.

\noindent $\bullet$ 用于给节点排序的判断函数: 
\verb|bool xUpsort(struct point p1, struct point p2)|.

\noindent $\bullet$ 用于打印向量元素的函数: 
\verb|void PrintVec(std::vector<double> vec)|.

% class Polynomial
\subsection{class Polynomial 的建构}
由于最终需要实现的程序需要完成生成插值多项式的任务，项目中不可避免地需要对多项式进行一个合理定义。本项目构建
了一个多项式类 \verb|class Polynomial|, 其具有多项式的一般性质，包括构造，求解，线性运算，并且能够打印输出。
该类存放于文件 \verb|Polynomial.h| 中，由于其是一个函数，首先我们对其继承\verb|Function|的虚函数以规范接口。
此外，其主要功能定义如下：

\noindent $\bullet$ 私有成员变量
\lstset{language=C}
\begin{lstlisting}
int deg; // 储存多项式次数
std::vector<double> coff; // 从低次到高次储存多项式系数
\end{lstlisting}
在此，我们把一个多项式的特征定义为次数(deg)和系数(coff)，通过对次数和系数的操作，我们就可以实现多项式
的所有操作。

\noindent $\bullet$ 构造函数
\begin{lstlisting}
Polynomial();   // 初始化构造，把多项式初始化为 0.
Polynomial(const std::vector<double> & vec);
                // 构造一个以 vec 为系数的多项式.
Polynomial(const Polynomial & poly);
                // 复制构造函数，通过深拷贝构造新的多项式类.
\end{lstlisting}

\noindent $\bullet$ 操作符重载
\begin{lstlisting}
double operator() (const double & _x) const;
        // 重载操作符()使其可以对输入值计算对应多项式值。
double operator[] (const int & n) const;
        // 重载操作符[]使其可以取出第n项系数，该操作是只读操作。
void operator= (Polynomial && poly);
        // 左值引用重载操作符 = 把给定 Polynomial 赋值给特定对象。
void operator= (Polynomial & poly);
        // 右值引用重载操作符 = 把给定 Polynomial 赋值给特定对象。
friend Polynomial operator+ (const Polynomial & p1,
                             const Polynomial & p2);
friend Polynomial operator+ (const Polynomial & poly,
                             const double & a);
friend Polynomial operator- (const Polynomial & p1,
                             const Polynomial & p2);
friend Polynomial operator- (const Polynomial & poly,
                             const double & a);
friend Polynomial operator* (const Polynomial & p1,
                             const Polynomial & p2);
friend Polynomial operator* (const double & c,
                             const Polynomial & poly);
// 友元函数重载操作符 +,-,* 使其完成普通多项式的加减乘法
// 以及对常数的运算
\end{lstlisting}
\textbf{附注1：} 在赋值运算 () 过程中，多项式乘法使用的是秦九韶算法，即
$$
a_0+a_1x+a_2x^2+\cdots+a_nx^n = a_0+x(a_1+x(a_2+x(\cdots(a_{n-1}+a_nx)\cdots)))
$$
\textbf{附注2：} 多项式乘法通过数组的乘法操作完成，其数学原理为
$$
\sum\limits_{i=0}^na_ix^i \sum\limits_{j=0}^mb_jx^j = 
\sum\limits_{k=0}^{n+m}c_k x^k, \quad c_k=\sum\limits_{i+j=k}a_i b_j
$$
\noindent $\bullet$ 功能函数
\begin{lstlisting}
double diff(const double & _x) const;
        // 覆盖继承多项式的求导函数运算
double second_diff(const double & _x) const;
        // 覆盖继承多项式的求二阶导函数运算
Polynomial GetDiff() const;
        // 把求导函数输出
int GetDeg() const; 
        // 获取次数
void assign(const std::vector<double> & vec);
        // 按照给定向量 vec 重构该多项式
void PrintCoff() const; // 系数形式输出多项式
void PrintPol() const;  // 多项式形式输出多项式
\end{lstlisting}
\textbf{附注:} 在求导运算过程中，多项式乘法使用的同样是是秦九韶算法。

\noindent 基于以上的建构，本项目实现了一个占据较小内存完成多项式功能模拟的函数类。

% class Bbasis
\subsection{class Bbasis}
类 \verb|class Bbasis| 是实现构造一个 B 样条基函数的类。其和根据节点构造的 B 样条基函数类都存在 \verb|Bbase.h| 这个
库函数中。函数的求值运算通过递归定义实现。

\bigskip
\noindent $\bullet$ 私有成员变量
\begin{verbatim}
int deg = 0;                       // 次数
int index = 0;                     // 所在基函数中的下标
const std::vector<double> * point; // 插值节点指针 下标超出部分   
\end{verbatim}
其中 \verb|point| 节点仅接收地址提高程序效率，下标用于定位在基函数中的下标，起始为0。

\bigskip
\noindent $\bullet$ 私有函数
\begin{verbatim}
double getpoint(const int & i) const // 根据下标返回插值节点的值
\end{verbatim}
下标在节点范围内按照节点输出，不在节点范围内则按照每次间隔 1 往两端取点。

\noindent $\bullet$ 构造函数
\begin{lstlisting}
Bbasis();                                    // 初始化构造
Bbasis(const int & _deg, const int & _index, // 构造函数
        const std::vector<double> & _point); // index 从 0 算起
\end{lstlisting}
\noindent $\bullet$ 公有成员函数
\begin{lstlisting}
double operator() (const double & _x); const // 对x求解插值函数值
double diff(const double & _x) const;        // 对x求解导数值
double second_diff(const double & _x); const // 求二阶导运算

int Getindex() const;  // 获得该函数的下标
\end{lstlisting}
\textbf{附注1：} 在赋值运算 () 和求导过程中，使用的都是递归定义，即
\begin{align*}
    B_i^{n+1}(x) & = \frac{x - t{i-1}}{t_{i+n} - t_{i-1}} B_i^n(x) +
    \frac{t_{i+n+1} - x}{t_{i+n+1} - t_i} B_{i+1}^n(x), \\
    \frac{d}{dx} B_i^n(x) & = \frac{nB_i^{n-1}x}{t_{i+n-1} - t_{i-1}}
    - \frac{nB_{i+1}^{n-1}x}{t_{i+n} - t_{i}}.
\end{align*}

% class Bbase
\subsection{class Bbase}

这个库旨在通过单个 B 样条函数 Bbasis 和一列插值节点构造对应于插值节点的 B 样条基函数。此外，由于 B 样条插值在实际
求解问题的时候用到的不止有 1-N 这 N 个基函数，因此此处允许对下标进行扩展，扩展的节点按照间隔1往节点两边延申。即 
$B_{-1}$ 对应的就是 $x_0-1$ 处构造的插值节点 B 样条函数。即设计信息如下：

\noindent $\bullet$ 私有成员变量
\begin{verbatim}
std::vector<double> point; // 插值节点
int deg = 0;               // 次数
std::vector<Bbasis> basis; // 生成的基函数 
\end{verbatim}

\noindent $\bullet$ 私有函数
\begin{verbatim}
void ComputeBasis() // 根据节点和次数生成基函数
\end{verbatim}

\noindent $\bullet$ 构造函数
\begin{lstlisting}
Bbase();                                   // 初始化构造
Bbase(const int & _deg,                    // 构造函数
      const std::vector<double> & _point); 
\end{lstlisting}
\noindent $\bullet$ 公有成员函数
\begin{lstlisting}
void operator= (Bbase && B)  // 重载运算符 "=" 完成赋值操作
void operator= (Bbase & B)
double operator() (const int & index,        // 根据下标指定基函数
                   const double & _x); const // 对x求解插值函数值
double diff(const int & index,
                   const double & _x) const; // 对x求解导数值
double second_diff(const int & index,
                   const double & _x); const // 求二阶导运算

int size() const   // 获取对应节点和次数所需基函数的个数
int GetDeg() cosnt // 获取次数
\end{lstlisting}

% class ppForm
\subsection{class ppForm}

ppForm 是一个通过多项式插值来构造对应节点多项式插值曲线的库，可以实现 $S_1^0$ 的线性插值类型和 $S_3^2$ 中
的五种线性边界条件。总体设计如下：

\noindent $\bullet$ 私有成员变量，用于唯一确定插值函数的信息
\begin{verbatim}
int type = S32;                   // 使用的插值类型，默认为 S32
std::vector<struct point> points; // 存储节点信息
std::vector<Polynomial> polyvec;  // 存储分段的多项式插值信息
// S32 的边值条件
int boundary = Nature;  // 边值条件，默认 Nature
double left = 0;        // 存储 x1 处初值条件
double right = 0;       // 存储 xN 处初值条件
\end{verbatim}

\noindent $\bullet$ 私有函数，用于辅助完成其他操作
\begin{verbatim}
void TypeForm();        // 根据给定参数进行求解器构造的函数
bool judgeXray() const; // 判断输入数据格式是否正确
int search_x(const double & _x) const; // 寻找 x 所在的区间
\end{verbatim}

\noindent $\bullet$ 构造函数
\begin{lstlisting}
ppForm(); // 初始化构造
// 按照节点信息构造
ppForm(const int & _type, 
       const std::vector<struct point> & _points, 
       const int & _boundary = Nature, 
       const double & _left = 0, const double & _right = 0);
// 按照节点向量信息构造
ppForm(const int & _type, const std::vector<double> & x, 
       const std::vector<double> & y,
       const int & _boundary = Nature, 
       const double & _left = 0, const double & _right = 0);
// 按照自变量值以及给定拟合函数进行构造
ppForm(const int & _type,
       const std::vector<double> & x, const Function & func,
       const int & _boundary = Nature);
\end{lstlisting}
特别的，对于给定待拟合函数的构造方式中，边界条件由 \verb|func| 直接计算得出。

\noindent $\bullet$ 公有成员函数，用于实现具体的插值函数计算功能
\begin{lstlisting}
double operator() (const double & _x); const // 对x求解插值函数值
double diff(const double & _x) const;        // 对x求解导数值
double second_diff(const double & _x); const // 求二阶导运算
void PrintPoly() const;        // 按照多项式形式打印存储多项式信息
void PrintPolyCoff() const;    // 按照系数形式打印存储多项式信息
void PrintPoint() const;       // 打印节点信息
\end{lstlisting}

\noindent $\bullet$ 确定每一段插值多项式的函数
\begin{lstlisting}
void ppFormS01();                 // 线性插值
void ppFormS32Complete();         // S32 Complete 类型
void ppFormS32Specified_Nature(); // S32 Specified 和 Nature 类型
void ppFormS32NotAKnot();         // S32 Not-a-kont 类型
void ppFormS32Periodic();         // S32 Periodic 类型
\end{lstlisting}
具体分析每一个函数的实现方式如下：

\bigskip
\noindent \textbf{\Rmnum{1}. S01 线性插值.}

只需要构建每一段的线性函数，利用关系
$$
s(x) - y_i = \frac{y_{i+1} - y_i}{x_{i+1} - x_i}(x - x_i)
$$
可以得到
$$
s(x) = \frac{y_{i+1} - y_i}{x_{i+1} - x_i} x + y_i -\frac{y_{i+1} - y_i}{x_{i+1} - x_i} x_i.
$$

\bigskip
求解 $S_3^2$ 类型前，首先明确总体流程。注意到对于三次插值样条，在 $[x_i,x_{i+1}]$ 这一段上可以表示为
$$
s(x) = f_i + s'(x_i)(x - x_i) + \frac{s''(x_i)}{2} (x - x_i)^2 + \frac{s'''(x_i)}{6} (x-x_i)^3.
$$
其中若记 $M_i = s''(x_i) \forall 1\leq i \leq N,$ 则具有递推式
\begin{align*}
    s'(x_i) & = f[x_i,x_{i+1}] - \frac{1}{6} (M_{i+1} + 2M_i)(x_{i+1} - x_i), \\
    s'''(x_i) & = \frac{M_{i+1} - M_i}{x_{i+1} - x_i}.
\end{align*}
同时对于二阶导具有以下关系：$\forall i = 2,3,\cdots,N-1$
$$
\mu_iM_{i-1} + 2M_i + \lambda_iM_{i+1} = 6f[x_{i-1},x_i,x_{i+1}],
$$
求解不同边值条件的 $S_3^2$ 类型插值，都可以按照一下流程进行：

% Step 1.
\noindent $\bullet$ 首先计算函数 $f(x)$ 对应插值节点的差分值 $f[x_i,x_{i+1}],f[x_{i-1},x_i,x_{i+1}].$ 
同时根据节点信息得到所有的 $\lambda_i, \mu_i$ 节点信息。

% Step 2.
\noindent $\bullet$ 根据边值条件和 $2-N-1$ 处二阶导关系可以建立线性方程组 $Ax = b,$ 求解得到的向量 
$x$ 存储的即为所有节点处二阶导信息。

% Step 3.
\noindent $\bullet$ 根据二阶导与一阶导，三阶导关系，应用对应关系式得到在节点 $1-N-1$ 处拟合曲线从右端
逼近的一阶导和三阶导数信息。

% Step 4.
\noindent $\bullet$ 根据导数和函数值信息构建唯一确定的多项式。

基于以上步骤，接下来不同边值条件的差别仅仅为求解关于二阶导数值的线性方程组 $Ax=b$ 不同而已了，以下仅仅
讨论需要的不同线性方程组。

\bigskip
\noindent \textbf{\Rmnum{2}. S32 Complete 类型.}

根据讲义 Lemma 3.6, 除了基本条件，还有条件
\begin{align*}
    2M_1+M_2 & = 6f[x_1,x_1,x_2], \\
    M_{N-1} + 2M_N & = 6f[x_{N-1}, x_N, x_{N}].
\end{align*}
由此可得
\begin{equation*}
    A = \left[
    \begin{matrix}
        2 & 1 & & & & & \\
        \mu_2 & 2 & \lambda_2 & & & & \\
        & & \ddots & & & & \\
        & & \mu_i & 2 & \lambda_i & & \\
        & & & & \ddots & & \\
        & & & & \mu_{N-1} & 2 & \lambda_{N-1} \\
        & & & & & 1 & 2
    \end{matrix}
    \right]
    b = \left[
    \begin{matrix}
        6f[x_1, x_1, x_2] \\ 6f[x_1,x_2,x_3] \\ \vdots \\
        6f[x_{i-1},x_i,x_{i+1}] \\ \vdots \\ 6f[x_{N-2},x_{N-1},x_N] \\
        6f[x_{N-1},x_N,x_N]
    \end{matrix}
    \right]
\end{equation*}

\noindent \textbf{\Rmnum{3}. S32 Specified 和 Nature 类型.}

其实这两种边值条件核心都在于给定了 $M_1$ 和 $M_{N}$ 两处的值，那么剩下 $N-2$ 个节点的二阶导函数值
正好可以由 $N-2$ 个原始的方程求得，则对应方程组为一个 $N-2$ 阶的方程组：
\begin{equation*}
A = \left[
\begin{matrix}
        2 & \lambda_2 & & & & & \\
        \mu_3 & 2 & \lambda_3 & & & & \\
        & & \ddots & & & & \\
        & & \mu_i & 2 & \lambda_i & & \\
        & & & & \ddots & & \\
        & & & & \mu_{N-2} & 2 & \lambda_{N-2} \\
        & & & & & \mu_{N-1} & 2
\end{matrix}
\right]
b = \left[
\begin{matrix}
        6f[x_1, x_2, x_3] - \mu_2 s''(a) \\ 6f[x_2,x_3,x_4] \\ \vdots \\
        6f[x_{i-1},x_i,x_{i+1}] \\ \vdots \\ 6f[x_{N-3},x_{N-2},x_{N-1}] \\
        6f[x_{N-2},x_{N-1},x_N] - \lambda_{N-1} s''(b)
\end{matrix}
\right]
\end{equation*}

\bigskip
\noindent \textbf{\Rmnum{4}. S32 Not a knot 类型.}

此处要求在 $x_2,x_{N-1}$ 处存在三阶导，即三阶导在此处连续即可。以 $x_2$ 为例，左右端
的三阶导分别满足
$$
    s'''(x_2)_- = \frac{M_1 - M_2}{x_1 - x_2}, \quad
    s'''(x_2)_+ = \frac{M_3 - M_2}{x_3 - x_2}.
$$
要三阶导存在，即要求 $s'''(x_2)_- = s'''(x_2)_+$, 计算可得
$$
(x_3 - x_2)M_1 - (x_3 - x_1) M_2 + (x_2 - x_1) M_3 = 0.
$$
同理有 $(x_N - x_{N-1})M_{N-2} - (x_{N} - x_{N-2}) M_{N-1} + (x_{N-1} - x_{N-2}) M_N = 0.$
对应方程组即为
\begin{equation*}
A = \left[
\begin{matrix}
        x_3-x_2 & -x_3+x_1 & x_2-x_1 & & & & \\
        \mu_2 & 2 & \lambda_2 & & & & \\
        & & \ddots & & & & \\
        & & \mu_i & 2 & \lambda_i & & \\
        & & & & \ddots & & \\
        & & & & \mu_{N-1} & 2 & \lambda_{N-1} \\
        & & & & x_{N}-x_{N-1} & -x_{N}+x_{N-2} & x_{N-1}-x_{N-2}
\end{matrix}
\right]
\end{equation*}
\begin{equation*}
b = \left[
\begin{matrix}
        0 \\ 6f[x_1,x_2,x_3] \\ \vdots \\
        6f[x_{i-1},x_i,x_{i+1}] \\ \vdots \\ 6f[x_{N-2},x_{N-1},x_N] \\
        0
\end{matrix}
\right]
\end{equation*}

\bigskip
\noindent \textbf{\Rmnum{5}. S32 Periodic 类型.}

周期类型在输入节点时要先判断首尾节点值是否相同。考虑可能由周期函数生成，而在计算机中不能算到绝对相同，
此处允许给定一个接近机器精度的误差范围。对于补充条件，即为补充两端点初一阶导，二阶导相同。即根据节点处
一阶导和二阶导的递推关系可得：

\begin{equation*}
A = \left[
\begin{matrix}
        1 & & & & & & -1 \\
        \mu_2 & 2 & \lambda_2 & & & & \\
        & & \ddots & & & & \\
        & & \mu_i & 2 & \lambda_i & & \\
        & & & & \ddots & & \\
        & & & & \mu_{N-1} & 2 & \lambda_{N-1} \\
        2(x_1 - x_0) & x_1 - x_0 & & & & x_N - x_{N-1} & 2(x_N - x_{N-1})
\end{matrix}
\right]
\end{equation*}
\begin{equation*}
b = \left[
\begin{matrix}
        0 \\ 6f[x_1,x_2,x_3] \\ \vdots \\
        6f[x_{i-1},x_i,x_{i+1}] \\ \vdots \\ 6f[x_{N-2},x_{N-1},x_N] \\
        6(f[x_1,x_2] - f[x_{N-1},x_N])
\end{matrix}
\right]
\end{equation*}
在实际求解时，为了防止出现极小的数据影响求解精度，求解过程中会对每行进行一定比例放大来提高求解精度。
此即为 ppForm 多项式插值类的设计介绍和原理。

% class Bspline
\subsection{class Bspline}

Bspline 类用于通过 B 样条插值进行计算，功能基本和 ppForm 相同，不同点在于本项目中该类未能实现 $S_3^2$ 
中的 Not a knot 的边值条件，而实现了 Cardinal 条件下的 $S_3^2,S_2^1$ 两种情况。其基本设计如下：

\bigskip
\noindent $\bullet$ 私有成员变量，用于唯一确定插值函数的信息
\begin{verbatim}
int type = S32;                   // 使用的插值类型，默认为 S32
std::vector<struct point> points; // 存储节点信息
Bbase basis;                      // 存储 B 样条基函数
std::vector<double> coeff;        // 存储系数
// S32 的边值条件
int boundary = Nature;  // 边值条件，默认 Nature
double left = 0;        // 存储 x1 处初值条件
double right = 0;       // 存储 xN 处初值条件
\end{verbatim}

\noindent $\bullet$ 私有函数，用于辅助完成其他操作
\begin{verbatim}
double getMul(const double & _a, const double & _b, 
              const double & _c) const; // 获取乘数
void TypeForm();        // 根据给定参数进行求解器构造的函数
bool judgeXray() const; // 判断输入数据格式是否正确
int search_x(const double & _x) const; // 寻找 x 所在的区间
\end{verbatim}

\noindent $\bullet$ 构造函数
\begin{lstlisting}
Bspline(); // 初始化构造
// 按照节点信息构造
Bspline(const int & _type, 
        const std::vector<struct point> & _points, 
        const int & _boundary = Nature, 
        const double & _left = 0, const double & _right = 0);
// 按照节点向量信息构造
Bspline(const int & _type, const std::vector<double> & x, 
        const std::vector<double> & y,
        const int & _boundary = Nature, 
        const double & _left = 0, const double & _right = 0);
// 按照自变量值以及给定拟合函数进行构造
Bspline(const int & _type,
        const std::vector<double> & x, const Function & func,
        const int & _boundary = Nature);
\end{lstlisting}
特别的，对于给定待拟合函数的构造方式中，边界条件由 \verb|func| 直接计算得出。

\bigskip
\noindent $\bullet$ 公有成员函数，用于实现具体的插值函数计算功能
\begin{lstlisting}
void operator= (Bspline && B); // 重载操作符 "="
void operator= (Bspline & B);
double operator() (const double & _x); const // 对x求解插值函数值
double diff(const double & _x) const;        // 对x求解导数值
double second_diff(const double & _x); const // 求二阶导运算
void PrintCoff() const;        // 打印系数
void PrintPoint() const;       // 打印节点信息
\end{lstlisting}

\noindent $\bullet$ 确定每一个 B 样条基函数对应系数的函数
\begin{lstlisting}
void BsplineS01();                 // 线性插值
void BsplineS32Complete();         // S32 Complete 类型
void BsplineS32Specified_Nature(); // S32 Specified 和 Nature 类型
void BsplineS32Periodic();         // S32 Periodic 类型

void BsplineCardinalS21();        // Cardinal S21 类型插值
void BsplineCardinalS32();        // Cardinal S32 类型插值
\end{lstlisting}
分析每一个函数的实现方式如下：

\noindent \textbf{\Rmnum{1}. S01 线性插值.}

根据基线性组合的公式，可得 $a_i = f_i$.

\noindent \textbf{\Rmnum{2}. S32 各种 类型.}

对于 $S_3^2$ 类型，根据基的线性组合，可以得到条件关系 $\forall i = 1, 2, \cdots, N,$
$$
f_i = a_{i-2}B_{i-2}^3(x_i) + a_{i-1}B_{i-1}^3(x_i) + a_{i}B_{i}^3(x_i).
$$
这里一共 N 个方程，而基函数有 N+2 个，包括 $B_{-1}^3,B_{0}^3,B_{1}^3,\cdot,B_{N}^3$. 
因此为了建立线性方程组 $Ax = b$, 还缺少 2 个条件。每种边值条件对应提供两个等式，按照基线性组合展开
即获得剩下的条件了。

\noindent \textbf{\Rmnum{3}. Cardinal 类型.}

这两种情况分别由 Theorem 47 (针对 $S_3^2$) 和 Theorem 48 (针对 $S_2^1$) 讨论完毕。注意 
$S_2^1$ 情况中计算函数值使用的是 $f(i+\frac{1}{2})$ 的值，而插值节点使用的仍然是整数节点。

% 项目测试
\section{项目测试}

具体应用前，使用两个库分别进行测试，主要测试其计算插值曲线的正确性以及对输入错误的检测能力。
为此，对于 $S_1^0$ 线性插值，以及 $S_3^2$ 三次插值中的 Complete, Specified, NotAKnot 三种
条件，都采用 Example 3.8 中的对数函数插值算例进行测试。而 $S_3^2$ 中的 Nature, Periodic 边值条件
则可以采用三角函数 $\sin x$ 来进行插值测试。

\subsection{ppForm 测试}

首先检测输入情况，结果存储到 testppResult.txt 中，给定节点 $(1,2), (-1,9), (0,4), (3,2)$ 进行
求解器生成，可以输出生成的的节点列为 $(-1,9), (0,4), (1,2), (3,2).$ 即已经根据节点自变量值完成了
排序，便于插值运算管理。

如果给定节点中自变量值出现重复，如给定 $(1,2), (-1,9), (0,4), (3,2), (1,3),$ 则会报错 
Values of x are wrong. 证明其可以对输入信息进行合理的判断。

此外，在计算时，插值函数只能计算给定插值节点之间的点的拟合函数值，对于区间之外的点计算时会报错。
如仍使用节点列 $(-1,9), (0,4), (1,2), (3,2),$ 使用线性插值测试，测试 $s(2)$ 可以正常输出插值的值 2, 
而测试 $s(4)$ 则报错 Out of range, 此时返回值默认为0.

接下来测试插值函数计算结果。把各种插值类型和边值条件分别使用对数函数和三角函数进行拟合得到拟合图像如下:

\begin{figure}[H]
\centering
\subfigure[$S_1^0.$]
{\includegraphics[width=4.8cm]{./figure/ppForm01.png}}          \label{pps01}
\subfigure[$S_3^2$ Complete.]
{\includegraphics[width=4.8cm]{./figure/ppForm32Complete.png}}  \label{ppComplete}
\subfigure[$S_3^2$ Specified.]
{\includegraphics[width=4.8cm]{./figure/ppForm32Specified.png}} \label{ppSpecified}
\subfigure[$S_3^2$ Nature.]
{\includegraphics[width=4.8cm]{./figure/ppForm32Nature.png}}    \label{ppNature}
\subfigure[$S_3^2$ Not a Knot.]
{\includegraphics[width=4.8cm]{./figure/ppForm32NotAKnot.png}}  \label{ppNotaKnot}
\subfigure[$S_3^2$ Periodic.]
{\includegraphics[width=4.8cm]{./figure/ppForm32Periodic.png}}  \label{ppPeriodic}
\caption{Test for ppForm} \label{Fig::ppForm}
\end{figure}

根据图像信息可以发现， ppForm 对于几种类型的插值都具有较好的拟合结果，可以认为其完成了插值拟合的要求。

\subsection{Bspline 的测试}

对 Bspline 类进行同样的测试，首先检测输入情况，结果存储到 testBsplineResult.txt 中，
给定节点 $(1,2), (-1,9), (0,4), (3,2)$ 计算，节点同样会被排序为 
$(-1,9), (0,4), (1,2), (3,2).$ 便于插值运算管理。

如果给定节点中自变量值出现重复，如给定 $(1,2), (-1,9), (0,4), (3,2), (1,3),$ 同样地会报错 
Values of x are wrong. 证明其可以对输入信息进行合理的判断。

此外, Bspline 对求解的函数定义域依旧有限制，对于插值节点区间之外的点计算时同样会报错。
如仍使用节点列 $(-1,9), (0,4), (1,2), (3,2),$ 使用线性插值测试，测试 $s(2)$ 可以正常输出插值的值 2, 
而测试 $s(4)$ 则报错 Out of range, 此时返回值默认为0.

接下来测试插值函数计算结果。把 Bspline 能实现的各种插值类型和边值条件分别使用对数函数和三角函数进行拟合
得到拟合图像如下:

\begin{figure}[H]
\centering
\subfigure[$S_1^0.$]
{\includegraphics[width=4.8cm]{./figure/Bspline01.png}}          \label{Bs01}
\subfigure[$S_3^2$ Complete.]
{\includegraphics[width=4.8cm]{./figure/Bspline32Complete.png}}  \label{BComplete}
\subfigure[$S_3^2$ Specified.]
{\includegraphics[width=4.8cm]{./figure/Bspline32Specified.png}} \label{BSpecified}
\subfigure[$S_3^2$ Nature.]
{\includegraphics[width=4.8cm]{./figure/Bspline32Nature.png}}    \label{BNature}
\subfigure[$S_3^2$ Periodic.]
{\includegraphics[width=4.8cm]{./figure/Bspline32Periodic.png}}  \label{BPeriodic}
\caption{Test for Bspline} \label{Fig::Bspline}
\end{figure}

同样根据图像可得 Bspline 针对经典的边值插值条件具有较好的拟合性质。针对 Cardinal 类型的插值测试在
后续解题中完成。

\section{项目求解题目应用}

\subsubsection*{A. Analyse the cubic-spline Interpolation of the function $f(x) = \frac{1}{1+25x^2}$.}

分别应用 ppForm 和 Bspline 对给定数量的点做等间距节点三次插值，绘图结果如下：

\begin{figure}[H]
\centering
\subfigure[ppForm]
{\includegraphics[width=7.2cm]{./figure/proAppForm.png}}   \label{ProAppForm}
\subfigure[Bspline]
{\includegraphics[width=7.2cm]{./figure/proABspline.png}}  \label{ProABsline}
\caption{Interpolate $f(x)$ by different knots.} \label{Fig::ProA}
\end{figure}

从图像可得，随着插值节点数量增多，曲线对 $f(x)$ 的拟合形态会变得越来越好。同时根据图像也可以发现这样的
插值有效地避免了 Runge 现象的出现，随着节点数增多，函数曲线整体性态都是拟合越来越好的。

进一步，以区间中点函数插值的最大值作为最大误差估计，
使用两种插值方式拟合的最大误差如下表所示：(原始数据在 result/ProAmaxerr.txt 内)

\begin{center}
\begin{tabular}{ccc}
    \toprule[1.5pt]
    \makebox[0.2\textwidth][c]{$N$} & \makebox[0.35\textwidth][c]{Maxerror of ppForm} 
    & \makebox[0.35\textwidth][c]{Maxerror of Bspline} \\
    \midrule[0.75pt]
    6  & 0.421705    & 0.411705    \\
    11 & 0.0205289   & 0.0205289   \\
    21 & 0.00316894  & 0.00316894  \\
    41 & 0.000275356 & 0.000275356 \\
    81 & 1.609e-005  & 1.609e-005  \\
    \bottomrule[1.5pt]
\end{tabular}
\end{center}

一方面，可以发现随之插值节点增多，误差基本上以每增长为 2 倍即可缩小至 $\frac{1}{10}$ 的速率进行收敛。
即误差收敛速率为幂函数阶，且两种插值方式收敛速度相同。

具体分析，设误差函数 $y=f(x),$ 由数据，函数大致满足关系
$f(2x) \approx 0.1f(x),$ 两边取对数可得 $\ln f(2x) = \ln(0.1) + \ln f(x),$ 一般的，可有 
$\ln f(2^t x) \approx t\ln(0.1) + \ln f(x) = \ln((0.1)^t f(x)).$ 取定 $x=1, s = 2^t,$ 
有 $f(s) \approx (0.1)^{\log_2 s}f(1) = s^{\log_2 0.1}f(1).$ 即最大误差收敛速度约为 
$O(N^{\log_2 0.1}) \approx O(N^{-3.32}).$
由此可以证明三次插值可以有效规避 Runge 现象出现。

\subsection*{B. Interpolate $f$ by the quadratic and cubic caldinal B-splines.}

实现方式见设计文档，调用的宏为 $S21Car(S_2^1 \mbox{ Cardinal}), S32Car(S_3^2 \mbox{ Cardinal}).$

\subsection*{C. Run the subroutines on the function $f(x) = \frac{1}{1+x^2}$.}

分别运用 Cardinal $S_2^1$ 和 Cardinal $S_3^2$ 对应的求解器运用给点节点对 $f(x),$ $x\in[-5,5]$, 
可得两种插值方式都能在规避 Runge 现象情况下较好地拟合曲线。绘图结果如下：

\begin{figure}[H]
\centering
\includegraphics[width=13cm]{./figure/ProC.png}
\caption{Use quadratic and cubic cardinal B-splines to interpolate.}
\label{Fig::Cardinal}
\end{figure}

\subsection*{D. Analyse the error of the interpolation.}

针对题目中给定的节点，带入 C 中的函数分别计算其和 Cardinal $S_2^1$ 和 Cardinal $S_3^2$ 
的差值，运算结果如下：

\begin{center}
\begin{tabular}{ccc}
    \toprule[1.5pt]
    \makebox[0.2\textwidth][c]{$x$} & \makebox[0.35\textwidth][c]{$E_s(x)$ of $S_2^1$} 
    & \makebox[0.35\textwidth][c]{$E_s(x)$ of $S_3^2$} \\
    \midrule[0.75pt]
    -3.5 & 6.93889e-017 & 0.000669568  \\
    -3   & 0.00141838   & 1.38778e-017 \\
    -0.5 & 1.11022e-016 & 0.0205289    \\
    0    & 0.120238     & 1.11022e-016 \\
    0.5  & 1.11022e-016 & 0.0205289    \\
    3    & 0.00141838   & 4.16334e-017 \\
    3.5  & 2.77556e-017 & 0.000669568  \\
    \bottomrule[1.5pt]
\end{tabular}
\end{center}

表中有一些数据接近机器精度，原因是这些点的函数值是对应情况的插值节点，理论上插值节点
处插值函数计算值应该等于拟合函数值，由于机器误差存在，会出现一个机器精度的误差偏移，
因此出现接近机器精度的数据。

更进一步，考虑非插值节点处的最大误差并结合图像，可以发现 $S_3^2$ 的拟合精度要比 $S_2^1$ 
更高。这里一方面是由于三次函数具有更高的光滑形态，拟合上可以适应更高阶导数；另一方面是
由于待拟合函数为偶函数， $S_3^2$ 在拟合时选取了 0 点，该点对于原函数也是一个形态比较突出
的点，在该处放置插值节点的拟合精度会比不放置插值节点的拟合程度更好。


\subsection*{E. Plot the heart.}

构建库函数 \verb|Heart.h| 用于计算拟合心形曲线，其求解思路总体是让曲线参数化。由于心形曲线方程满足
$$
x^2 + (\frac{3}{2}y - \sqrt{|x|})^2 = 3.
$$
由于图像左右对称，可以做参数化，对 $\forall \theta \in [0, 2\pi),$
\begin{align*}
    x & = \sqrt{3} \sin \theta \\
    \frac{3}{2} y - \sqrt{|x|} & = \sqrt{3} \cos \theta
\end{align*}
经过化简得到
\begin{align*}
    x & = \sqrt{3} \sin \theta \\
    y & = \frac{2}{3} (\sqrt{\sqrt{3}|\sin \theta|} + \sqrt{3}\cos \theta)
\end{align*}
做伸缩变化 $\theta = 2\pi t,$ 则参数化曲线 $\gamma(t) = (x(t), y(t)):[0,1)\to \mathbb{R}^2$ 
即为心形曲线的参数方程。应用 cumulative chordal lengths 来获取各个参数值对曲线的拟合插值点，
可以通过构造两个 Bspline 求解器来分别插值拟合 $x,y$ 两个坐标进而获取总体插值曲线信息。

理论分析心形曲线，其有两个点 $(0,\frac{2}{3}\sqrt{3}),(0,-\frac{2}{3}\sqrt{3})$ 是尖点，此处导数
无定义。因此在选取插值节点时，选取其中一个作为起点和中点，另一个点作为插值节点中间节点，其他点按照左右
对称与 $y$ 轴等角度间距取点即可。考虑到两个尖点处的奇异性，这两个点附近的插值节点选得更加近一些。
边值条件方面，由于起点中点都无导数意义，在此统一选取 Nature 类型作为边值条件以减少束缚。
对 $N=10,40,160$ 分别进行插值计算，绘图结果如下：

\begin{figure}[H]
    \includegraphics[width=13cm]{./figure/ProE.png}
    \caption{Interpolation to heart curve by different number of knots.}
    \label{Fig::heart}
\end{figure}

从图像可以发现，随着插值节点增多，插值曲线对心形曲线的拟合形态会逐渐变好，对于奇异的两个点的拟合
也会随着节点数增多而逼近。

\section{总结}

该项目实现了 ppForm 和 B-splines 插值的多种方式，可以用于进行一些简单的插值函数拟合，
具有数值逼近的意义。

\end{document}